\chapter{非线性方程的数值解法}
\section{单方程}
如果函数$ f $是连续的，且$ f(a)<0<f(b) $，其中$ a<b $，那么在一定的条件下，零值定理告诉我们在$ (a,b) $上一定存在一个点使得函数$ f $为0。

此时可用\textbf{二分法}，即首先确定两个点，使得函数值一正一负，然后在这两个点之间进行二分。

二分法收敛比较慢，但它可以确保收敛。\textbf{牛顿方法}收敛速度比较快，但并不一直收敛。该方法的思路在于对于任何非线性函数$ f(x) $，我们用它的切线来近似$ f $，然后该切线的与X轴的交点可以作为$ f $的根。可以看到如果$ f $是线性函数，这种方法是一次迭代就搞定。如果$ f $是非线性的，则需要反复迭代。

依据该思路，对于任何初始猜测$ x_k $，它对应的切线上的点坐标为$ (x_k,f'(x_k)) $，该切线的方程就可以写为，
\[ y-f(x_k)=f'(x_k)(x-x_k) \]

该切线与X轴的交点的横坐标可以令$ y=0 $而得到，因此，有\footnote{也可以利用泰勒一阶展开来理解该公式。}，
\begin{equation}\label{ne_eq_1}
x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}	
\end{equation}

另外，还有一种\textbf{割线法}，在大规模方程求解中，可能会极大加快速度。它利用函数上的两点形成的割线的斜率作为\eqref{ne_eq_1}式中斜率的替代，即，
\[ x_{k+2}=x_k-\frac{x_{k+1}-x_k}{f(x_{k+1})-f(x_k)}f(x_k) \]
这里，$ x_{k+2} $就是要计算的点，$ x_k,x_{k+1} $是曲线上的两个点。可以证明，牛顿法和割线法最终都会收敛到真值。

\section{方程组}
\subsection{牛顿方法}
比如有如下方程组，
\begin{align*}
	f_1(u,v,w)=0\\
	f_2(u,v,w)=0\\
	f_3(u,v,w)=0
\end{align*}
可以定义向量值函数$ F(u,v,w)=(f_1,f_2,f_3)' $，则上述方程组可以简写为$ F(\bm{x})=0$，其中,$ \bm{x} = (u,v,w)' $。那么对应于单方程牛顿法的$ f' $，此时是一个雅可比矩阵，定义为，
\[ DF(\bm{x})=\begin{bmatrix}
	\frac{\partial f_1}{\partial u}& 	\frac{\partial f_1}{\partial v}& 	\frac{\partial f_1}{\partial w}\\
	\frac{\partial f_2}{\partial u}& 	\frac{\partial f_2}{\partial v}& 	\frac{\partial f_2}{\partial w}\\
	\frac{\partial f_3}{\partial u}& 	\frac{\partial f_3}{\partial v}& 	\frac{\partial f_3}{\partial w}\\
\end{bmatrix} \]
对应于\eqref{ne_eq_1}式的迭代公式此时为，
\[ \bm{x}_{k+1}=\bm{x}_k-(DF(\bm{x}_k))^{-1}F(\bm{x}_k) \]

考虑到矩阵求逆比较耗时，可以进一步使用求解线性方程组的高斯消去法来得到矩阵的逆。
\subsection{Broyden方法}
对于一个900个变量900个方程的系统，牛顿法每步需要计算$ 2\times 900^2+900 $个标量函数值，运算代价是高昂的。因此，可以对单方程中的割线做一个多变量推广，缓解对雅可比矩阵的计算压力。

假设$ \bm{A}_k\equiv \bm{J}(\bm{x}_k)  $是第$ k $步迭代估计的雅可比矩阵，$ \bm{w}_{k+1}=\bm{x}_{k+1}-\bm{x}_k $是步长，则函数的变化为$ \bm{y}_{k+1}=f(\bm{x}_{k+1})-f(\bm{x}_k) $。很明显，因变量之差等于切线斜率乘以自变量之差，即，
\begin{equation}\label{ne_eq_b1}
\bm{y}_{k+1}=\bm{A}_{k+1}\bm{w}_{k+1}	
\end{equation}

我们的目的在于得到$ \bm{A}_{k+1} $，而$ \bm{A}_{k+1} $是一个方阵，有$ n^2 $个元素。而上面的公式意味着只有$ n $个方程，尚不足以得到全部的$ \bm{A}_{k+1} $中的元素。Dennis and Schnabel (1983)表明可以再加上如下条件，
\begin{equation}\label{ne_eq_b2}
\bm{A}_{k+1}\bm{z}=\bm{A}_k\bm{z},\qquad \text{with}\qquad (\bm{x}_{k+1}-\bm{x}_k)'\bm{z}=0	
\end{equation}
其中$ \bm{z} $是任意向量。\eqref{ne_eq_b1}式和\eqref{ne_eq_b2}式保证了函数$ f(\bm{x}) $在点$ \bm{x}_k $和$ \bm{x}_{k+1} $的相邻线性近似的差是最小的。

\eqref{ne_eq_b2}式表明，
\[ (\bm{A}_{k+1}-\bm{A}_k)\bm{z}=\bm{0}, (\bm{x}_{k+1}-\bm{x}_k)'\bm{z}=0\]

这表明对于任意的向量$ \bm{u} $，下式总成立，
\begin{equation}\label{ne_eq_b3}
	\bm{A}_{k+1}-\bm{A}_k=\bm{u}(\bm{x}_{k+1}-\bm{x}_k)'
\end{equation}

含义就是$ \bm{A}_{k+1}-\bm{A}_k $中的每一行与$(\bm{x}_{k+1}-\bm{x}_k)'  $这一行成常数比例。那么接下来的推导就显而易见，
\[\begin{aligned}
(\bm{A}_{k+1}-\bm{A}_k)(\bm{x}_{k+1}-\bm{x}_k)&=\bm{u}(\bm{x}_{k+1}-\bm{x}_k)'(\bm{x}_{k+1}-\bm{x}_k)\\
\Longrightarrow \bm{y}_{k+1}-\bm{A}_k(\bm{x}_{k+1}-\bm{x}_k)&=\bm{u}(\bm{x}_{k+1}-\bm{x}_k)'(\bm{x}_{k+1}-\bm{x}_k)\\
\Longrightarrow \bm{u}&=\frac{\bm{y}_{k+1}-\bm{A}_k(\bm{x}_{k+1}-\bm{x}_k)}{(\bm{x}_{k+1}-\bm{x}_k)'(\bm{x}_{k+1}-\bm{x}_k)}\\
\Longrightarrow \bm{A}_{k+1}&=\bm{A}_k+\frac{(\bm{y}_{k+1}-\bm{A}_k\bm{w}_{k+1})\bm{w}_{k+1}^T}{\bm{w}_{k+1}^T\bm{w}_{k+1}}\hspace{2em}\text{利用了\eqref{ne_eq_b3}式}
\end{aligned}
\]

从中可以看到，
\begin{itemize}
	\item 计算量大大减少了。
	\item 有时即便是初始的雅可比矩阵也计算困难，则可直接用单位矩阵替代。
\end{itemize}

